function [nlml,hyp_prior,hyp_post,ytest] = Run_icu24gp(x,y,xtest,DOPLOT)
%RUN_ICU24GP	compute GP analysis for the first 24 hours after ICU admission
%	[hyp_prior,hyp_post,ytest] = Run_icu24gp(xtrain,ytrain,xtest,DOPLOT) 
%	
%
%	Inputs:     x, y                    input, output variables for
%                                       training
%               xtest                   input test variable
%               DOPLOT                  1, represent plot
%		
%		
%
%	Outputs:    hyp_prior, hyp_post     hyperparameters structure prior and
%                                       posteriors
%               ytest                   double-column vector with output mean and
%                                       variance

%		
%		
%	See also Gui_icu24gp

%	References: N/A
%	
%	

%	Copyright 2013 MAF Pimentel
%	This program is free software: you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation, either version 3 of the License, or
%	(at your option) any later version.
%	
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%	
%	You should have received a copy of the GNU General Public License
%	along with this program.  If not, see <http://www.gnu.org/licenses/>.


%	$Author: MAF Pimentel$
%	$Revision: 1.0.0.0$
%	$Date: 12-Jun-2013 15:36:05$
%	Contact: marco.and.pimentel@gmail.com
%	Originally written on: PCWIN64

%% Apply gaussian processes to time series using NLML-based training
% Define covariance function
k1 = @covSEiso;                         % long term trend
k2 = @covSEiso;                         % short term length-scales
covfunc = {@covSum, {k1, k2}};          % sum up the covariance terms

hyp.cov = [3 0 1 0]; hyp.lik = -2;      % initialise hypers

% Apply timeseries normalisation
y1 = (1/std(y))*(y - mean(y));

% Define which hyperparameters to train
params = [3 1 3 .3; ...  % which parameter to optimise; lower-end of interval; upper-end of interval; step in interval;
          1 -2 1 .4];

% Train GP
hyp.cov = gptrain(hyp,covfunc,x,y1,params);     % optimise parameters
hyp_prior = hyp;                                % save hyp prior
[hyp, ~, ~] = ...
        minimize(hyp, @gp, -200, @infExact, [], covfunc, @likGauss, x, y1);
nlml = gp(hyp, @infExact, [], covfunc, @likGauss, x, y1);
hyp_post = hyp;                                 % save hyp posterior
[mu, s2] = gp(hyp, @infExact, [], covfunc, @likGauss, x, y1, xtest);

ytest = [(mu*std(y) + mean(y)) s2*std(y)];

if DOPLOT
    figure;
    zz = xtest; mu = ytest(:,1); s2 = ytest(:,2);
    f = [mu+2*sqrt(s2); flipdim(mu-2*sqrt(s2),1)];
    fill([zz; flipdim(zz,1)], f, [7 7 7]/8); alpha(0.25); hold on;
    plot(zz, mu,'--','Color','k'); plot(x,y,'*','Color','r');
    hold off; xlabel('x, input'); ylabel('y, output');
    title('GP-Regression');
end

end